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FS U INTRODUCTTON 


As sound waves propagate through the ocean, the complex 
acoustic pressure field which is generated by the source 
depends mainly on the path followed by the acoustic rays and 
the sound speed along this particular path. Due to this 
close relationship between the sound speed field and the 
acoustic pressure field, an attempt may be made to estimate 
the range-dependent sound speed profile (SSP) between a 
fixed source and an array of receivers. The characteriza- 
tion of the SSP from acoustical measurements generally 
involves inverse techniques in order to infer the acoustic 
properties of the medium from the pressure field measured at 
the receivers. 

Due to the complexity of the ocean the inverse problem 
is most often non-linear and underdetermined. Classicial 
acoustic tomography solves this problem by linearization. 
The tomographic method is able to estimate the perturbations 
of the sound speed field by comparing the measured travel 
times of particular rays with those computed numerically 
from a reference sound speed field and a raytrace or a 
normal mode algorithm (Munk and Wunsch, 1978). The 
procedure provides maps of the perturbations in the sound 
speed field and indicates how different the actual field is 


from the one used as a reference. If a large enough number 


of ray paths are used in the computations, the spatial 
resolution in the map of the sound speed field may be better 
than the one obtained from discrete CTD measurements (Howe, 
1986). 

Matched field processing is a different type of inverse 
method which was first proposed as a method to locate an 
acoustic source in the ocean. The principle is to compare 
the measured complex acoustic pressures at a vertical array 
with those computed from an acoustic model using various 
positions of the target source (Bucker, Я Тре 
procedure generates a function which is a measure of the 
likelihood between the actual acoustic pressure field 
created by the source (unknown position) and a replica 
pressure field generated from an estimate of the source 
location: 

This study is an attempt to use matched field processing 
as an alternate tool to solve the inverse problem in 
acoustic tomography. Given the position of the source and 
the receivers, matched field processing 15 employed to 
compare the true complex acoustic pressures at the receiving 
array with the ones computed from an acoustic propagation 
model and various sound speed fields. Computer simulations 
are used to demonstrate the performance of various 
estimators under different signal-to-noise ratios and noise 


correlation matrix Structures. 


Chapter II provides a theoretical presentation of the 
likelihood estimators which are used in this study. It also 
illustrates how the noise is modeled and added to the 
Simulated data. The simulations are shown in Chapter III 
for various conditions of noise. The first simulations deal 
with the simple case where only one parameter of the medium 
is unknown (e.g., SOFAR axis depth, surface sound speed, 
Single frontal boundary) and where the noise is absent. The 
next cases are applied to the localization of an eddy in a 
noisy medium. Situations of both spatially uncorrelated 
noise and correlated noise were examined. Comparison of the 
estimators is provided in Chapter IV; the spreading of each 
likelihood function about the true value is examined in more 
detail under several conditions of noise power and noise 
correlation. Also discussed is the problem of incomplete or 
poor knowledge of the other environmental parameters, e.g., 
incorrectly specifying the bottom absorption property, and 
how this may introduce inconsistency in the procedure and 
lead to a decrease in the efficiency of the likelihood 


то Стопе. 


II. ACOUSTIC TOMOGRAPHY ANDIMATCHED ETELDS 


А. PRINCIPLES OF MATCHED FIELD PROCESSING 

Classical beamforming for plane waves is obtained by 
measuring the maximum likelihood between the actual value of 
the complex signal at each hydrophone and the values 
computed from an expected bearing (Ziomek, 1985). 

In a similar fashion, the distance to a target in the 
near field can be estimated by comparing the actual values 
with those computed for different distances. The range is 
assumed to be correct when both sets of values correspond to 
the same wave front curvature (Ziomek, 1985). 

Matched field processing has been used traditionally to 
find the location of an acoustic source in a shallow water 
environment. The general principle is to store the values 
of the received signal (amplitudes and phase) at each 
element of the array and then compare them with theoretical 
values computed for different possible positions of the 
target source. The true location of the emitter is 
determined when both fields match (Bucker, 1976; Baggeroer, 
Kuperman and Schmidt, 1988). 

Different criteria may be used to measure the likelihood 
or degree of matching. Each one generates a different 


function which is generally well adapted for a particular 


type of noise. Matched field detection is consistent when 
the medium is completely determined. 

Matched field tomography deals with the inverse problem. 
Given that the source location is known, the purpose of the 
procedure 15 Жо estimate the medium characteristics, 
particularly the range-dependent sound speed profile. Tn 
this case, the replica or estimated fields are built from 
many sound speed profiles and one tries to match them with 


the measured one. 


В. W r Py INSJNOISE-PREE CONDITIONS 
I Bucker Method 
According to Bucker (1976), the following "detection 
factor" may be used as a measure of the difference between 
the exact and the replica fields: 


NR NR 
«ко, ве KR. 
ee _ y: jk jk 

Buck = 2=1 14) 


Ñ 2s) 


where terms are defined as follows: 
KR = spatial autocorrelation matrix of one replica 
field, 


KS = spatial autocorrelation matrix of the actual 
fleld, 


NR = number of hydrophones in the array, 


F = scaling factor to insure a result between 
0 and 1, 


<> = time average used to remove the component due 
to the noise, 


= complex conjugate 


The matrices KS and KR are defined by 


KSjk 7 BjBk- (2.2) 


КВ-к 


P'jBP'x- (2.3) 


where pj and р'- denote, respectively, the complex envelope 
of the acoustic pressure and the replica pressure at 
hydrophone j; similarly for hydrophone k. As demonstrated 
by Bucker (1976), it is convenient to define the correlation 
matrices from the complex envelopes of the signal because 
the rapidly-varying time component is removed from the 
computations. These complex envelopes are easily obtained 
by processing the incoming signal through а classical 
quadrature demodulator (multiplication by a sine wave 
followed by a low-pass filter). 

In the absence of noise, the time average is 
unnecessary (KS is time independent) and the normalized 


detection factor becomes: 


NR МЕ 
) ) КО. „КК. 
жк J “~ _ ____ | 
век Жы Е 12 № МЕ 1/2 (2.4) 
(S у KR. KR UN y Кок Ка“) 
j=1 k=12j J J 451 k=Î#) 


реши се = Similar to the classical correlation 
coefficient of two random variables. The expression given 
above does not use the diagonal elements of the matrices and 
can be interpreted as the output of a regular beamformer. 
Its value is one in the case of complete equality of the 
fields (KR = KS). We note also that the Bucker detection 
factor is one when the matrices KR and KS are proportional, 
аз а consequence of the Schwarz inequality. 

2. Heitmeyer Method 

Heitmeyer (1984) defined another detection factor 

which he called the "source location ambiguity function." 


Its value is given by the following expression: 
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yey En En | 
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and may be rewritten in terms of phases and magnitudes: 
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It can easily be seen that the function is 


unnormalized. From the inequality: 
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The ambiguity function is always bounded by a 


quantity that may be considered as the average power 


detected at each hydrophone. 


When the replica and the measured fields match 


exactly, the expression reduces to the following: 


р' = Bn (2258 
which yields 
Fl 
HEIT = nel x 5 3 xd : (2-15) 
Ы NR 2 CONES = 4 
NR ) P. 
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The inequality demonstrated previously becomes an identity 


when the actual and replica fields are identical. 


In order to obtain a normalized ambiguity function, 


we divide Equation (2.6) by this upper limit. 


| NR | 2 
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3. Relation between Bucker and Heitmeyer Methods 


Using the definition of the spatial autocorrelation 


Matrix, Equation (2.4) may be written: 
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If we allow the subscripts j and k to be equal, this 


expression becomes: 


hence, 


КЕ NR 
= pi ¿ PkPk 

ВОСК = = ; NR 5 (2.14) 
1р." ) ру 
j=1 J =1 


or, by changing the subscripts, 


BUCK = (2815) 
у 


which is exactly the normalized ambiguity function defined 
by Hertmeyer іп Equation (2 11). 

More generally, Бу developing the expressions of 
both functions, we can derive the following relation between 


the Bucker and Heitmeyer definitions: 


NR NR 
с” v | 1-1 
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ыг се ищ 26 
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BUCK = HEIT > 


This relation is not a simple proportionality raci; 
because it changes with the replica fields p's. 
4. Center of Gravity Method 
Other likelihood estimations can be developed using 


any function which has a maximum in case of perfect 
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matching. The center of gravity method is a procedure which 
solves the problem from a more mathematical viewpoint. 

In the complex plane, vectors P and P' define two 
sets of NR points, where NR represents the number of 
receivers of the array. The components of the points are 
the real and imaginary parts of the complex acoustic 


pressures: 


— 


ЕЕ is composed Ci points (pj, cos jPi Sin +.) 


а с роте ир Cos") ,p"; sin +š) 


In order to compute a detection factor, we calculate 
the euclidean distance between the centers of gravity G and 
G' of both sets (Figure 2.1). This distance is inversely 
related to the likelihood of the fields. 

With a sum of weights of 1, the points G and G' are 


given by their coordinates: 


pes 
Xg = NE LN COS 133 (2227) 
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The distance becomes: 


D = ее” + Е р. (2. 2 
ОТ 
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We then normalize the quantity in order to obtain 


Unity for complete таван ој 


DN 35525 (2.29 


where Dgà4, is the largest unnormalized distance among all 
the replicas. 

DN is only an estimation of likelihood. Although it 
is possible for two different sets of points to have the 
same center of gravity, if they are concentric, this does 
not occur in the simulations and the method keeps its 
consistency. We will see later that this distance function 
may lead to high secondary lobes and thus is not always 


reliable: 


C. THEORY IN PRESENCE OF NOTSE 
Although it is possible to use the former expressions 


when noise is present which contaminates the signal, the 
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following two functions are more specifically suited for use 
in the presence of noise. 

Johnson (1982) previously demonstrated the equivalence 
between the problem of bearing determination and the 
estimation of the spectrum of a signal. Due to this 
Similarity, all modern spectral estimation algorithms apply 
equally to beamforming and matched field processing. 

Although many functions may be used, as for example, 
Meee (Schmidt, 1981) or linear predictor (Johnson, 1982), 
we will particularly emphasize the Bartlett and Maximum 
Likelihood parameters which are two powerful estimators in 
eaget location problems. These estimators are especially 
useful in noisy conditions, because Baggeroer and his 
colleagues (1988) showed that they reduce to the Bucker or 
Heitmeyer structures when the noise is absent. 

1. Bartlett Method 

This method comes directly from spectral estimation 
theory. Baggeroer, Kuperman and Schmidt (1988) demonstrated 
that the power output of a Bartlett beamformer could be 


written in the following quadratic form: 


BART 


W KTW (2224) 
where W represents the normalized velocity potential vector 
of the replica field and KT is the total spatial correlation 


matrix of the signal embedded in noise. 
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Due to the proportionality between the velocity 
potential and the pressure field, the equivalent expression 


will be used: 
BART - (P'/|P'|)* KT (P'/|P']) (2.25) 


where р! is the complex acoustic pressure vector of the 
replica at the array. Under the condition of perfect 
matching and no noise, 

ое M x T 

BART? =.) Von) sae (2.26) 

1=1 7 
which is the summation of all signal powers among the 
hydrophones. 

In order to normalize the function, we will divide 
the estimator by its largest value. For comparison between 
the different methods, we will focus on the width of the 
main lobe rather than its absolute value. 

ТЕ рі апа пу епос respectively, the signal and 
the noise pressure at hydrophone i, the spatial correlation 


matrix has the value: 


КТ а = Е((рү+п{) (ру+п-)”) (2.25) 


where E( ) denotes the expectation operation. 
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When the signal and the noise are uncorrelated, the 
matrix reduces to the simple sum of signal and noise 


matrices: 


КТіз - E(PiP$3') + Е(пјпј') (2.28) 


KTij = KS1j + KN i 5 (25:29) 


These matrices are hermitian and at least 
semidefinite. 
2. Maximum Likelihood Method 
In spectral estimation, the Maximum Likelihood 
method, also called Capon's method or the minimum variance 


method, is used to compte the power spectral density of a 


random process (Kay, 1988). Its expression is given by: 

Puy (f) = (eH R," 1 e)71 (2.30) 
where: 

f = frequency, 

Ryx = time correlation matrix of the process, 

e = vector whose it component is eJ? f, 

H = transposition of the conjugate matrix. 


In a similar way the output of a Maximum Likelihood 


beamformer is defined: 


1:56 


MIM = (W: кт 1 и) -1 (2.305 


where W and KT have previously been defined. 


Following the procedure of Equation (2.25), Equation 


(2.31) 1s modified: 
MLM = ((P'-/|P']) KT71 (P'7|P'|))71 201 


In the absence of noise, KT reduces to the spatial 
autocorrelation of signal only, KS (Equation (2 29 пи As 
can be seen when the array is composed of two hydrophones, 
the matrix is generally singular and has no inverse. The 
calculation is made possible by adding a small amount of 
noise to the diagonal. As with the Bartlett estimator, the 
Maximum Likelihood factor will be divided by its largest 


value for normalization. 
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Figure 2.1 Distance Between Two Complex Sets of Points 


ТТТ. SIMULATION OF ACOUSTIC TOMOGRAPHY 


Matched field processing has the same form as classical 
beamforming. However, instead of comparing the actual field 
vector with a plane wave replica vector, we will try to 
match the actual vector with a vector computed from an 
acoustic propagation model. Measured data will also Бе 
simulated with the same code, then embedded or not in noise 


depending on the scenario under investigation. 


А. PROCEDURE 
1. Description of the Simulation 

The receiver is modeled as a vertical array and is 
assumed to be composed of 20 hydrophones evenly distributed 
between 550 m and 1500 m. The source is located 100 km from 
the receiver at a depth of 1000 m. It emits a pure sine 
wave (tonal) centered at 100 Hz. Due to the inherent 
limitation of vertical angles in the parabolic approximation 
(Ziomek, 1985), the transmitter was selected to have a 
beamwidth of 40°. The bottom is 5000 m deep and is assumed 
to be flat and fully absorbing. This choice was made to 
speed up the calculations and is not a restrictive 
assumption. It assumes all the energy propagates by water- 


borne paths. 
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2. Acoustic Propagation Model 

Several models could have been used to simulate the 
acoustic fields, for example, normal mode theory or the 
parabolic equation (PE). The PE model was used because it 
1s more suitable for a deep water simulation; for example, 
an ocean bottom of 5000 m allows almost 670 propagating 
modes at 100 Hz and would have been computationally 
intensive using normal mode theory. 

The classical PE approximation with  split-step 
Fourier transform (Coppens, 1982) was available in the 
Environmental Acoustic Research Group package of models 
resident at NPS. The source code was slightly modified to 
save the complex acoustic pressures at the hydrophones in a 
file. The measured and the replica complex pressure fields 
were then stored in order to run the simulation programs. 

3. Simulation of Noise 

Noise was added to the measured data in order to 
produce a realistic problem and to study the behavior of the 
estimators in different environments. Following the 
procedure described by Porter, Dicus and Fizell (1987), 
noise was introduced by the mean of its spatial correlation 
matrix. This procedure is better than just altering the 
data with random noise because the probability density 
function of the noise is difficult to estimate. Moreover, 
all the estimators considered were written in terms of 


Correlation matrices. 
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Ambient noise falls in two categories: 

- uncorrelated noise. Its correlation matrix 15 
proportional to the identity matrix ала gene 
proportionality factor is an indicator of the noise 
power. 

- correlated noise. In this case the matrix has non zero 
terms outside the diagonal, but is nevertheless an 
hermitian matrix. 

Several attempts to measure the coherence of ambient 
noise in the ocean have been conducted during the past 
years. One of them was made from the Trident Vertical Array 
and is described by Urick (1984). Figure 3.1 depicts the 
results of this study and has been used to generate a model 


of the noise correlation matrix 


The matrix was modeled in the following way: 
KNij = 52 2) 0-21 (3.1) 


where :? depends on the noise power and > is a factor which 
indicates how fast the coherence falls off outside the 
diagonal. Тһе larger 8 is, the more uncorrelated the noise 
is, i.e., its spatial correlation scale becomes shorter. 

Due to the spacing between the receivers in the 
array and the frequency (100 Hz) used in the simulation, 
Figure 3.1 shows that £= 1.7 is consistent with the 
observed ambient noise correlation. 

Although the spatial correlation matrix of the 
noise, КМ, 15 generally a complex hermitian matrix, the 


analytic modeling shown in Equation (3.1) describes a real 
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symmetric matrix. As explained by Cox (1973), this 
approximation is valid in the special case of zero time 
delay, i.e., when it is assumed that the noise is in phase 
among all the hydrophones of the vertical array. This 
assumption is relatively consistent for low frequency noise. 
Below 150 Hz the noise iS principally due to distant 


shipping and arrives mainly from the horizontal. 


B. ONE DIMENSIONAL PROBLEM (NOISE-FREE CONDITIONS) 

In the following simulation, the shape of the sound 
speed profile is the only unknown. If the profile is 
digitized in 1 m intervals, then for a water depth of 5000 
m, one would have to determine 5000 values to match the 
complete sound speed profile. Such a procedure would lead 
to an unmanageable number of computations, especially if one 
tries to match a large number of replica fields. However, 
as we are only interested in demonstrating the feasibility 
of the procedure, we will begin by studying the simple cases 
where only one or two points of the sound speed profile are 
unknown. 

1. Determination of a SOFAR Axis Depth 

We initially start with a bi-gradient sound speed 
profile having a sound speed minimum at 1000 m (Figure 3.2). 
Two replica profiles are considered wherein the SOFAR axis 
depth is altered by +/- 200m (step 10 m). Note that 


because the surface and bottom sound speeds are unchanged, 


zu 


the change in axial depth results in a change in the 
gradients of both the upper and lower segments of the SSP. 

Three estimation techniques, the Bucker, Heitmeyer, 
and center of gravity methods are utilized to determine the 
true depth of the sound speed minimum. These estimators are 
computed from Equations (2.4), (2. 11) ana n The 
estimated depth of the SOFAR axis is found when an estimator 
shows a peak with a detection factor of 1. The results are 
shown in Figures 3.3, 3.4 and 3.5; each estimator 
demonstrates a different behavior. 

The Bucker detection factor and the Heitmeyer 
ambiguity function both indicate a maximum at the true 
location of 1000 m. However a strong side lobe, centered at 
1040 m, indicates these two detection factors are not robust 
enough to provide an unambiguous selection of the SOFAR axis 
depth. In addition, strong secondary side lobes are also 
present. 

The center of gravity method (Figure 3.5) proves to 
be a better estimator for this situation. The main lobe is 
much narrower and no other lobes exist. For a noise-free 
ocean, this is the best estimator among the three to solve 
this particular problem. 

2. Determination of Surface Sound Speed 

In order to mimic typical seasonal ог spatial 

changes in the SSP, the surface sound speed was permitted to 


fluctuate by */- 5 m/s about a mean value (step 1 m/s). For 
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this situation the SOFAR axial depth was fixed at 1000 m. 
Hence, only the upper gradient changes as seen in Figure 
3.6. Ме seek an estimator that will match the true surface 
sound speed (solid line in Figure 3.6). Using the same 
three estimation techniques as above, the results are shown 
in Figures 3.7, 3.8 and 3.9. For this situation, we observe 
a nearly identical behavior of the Bucker and the Heitmeyer 
functions, both of which have moderate side lobes at about 
0.7. As before, the center of gravity estimator remains the 
best without any ambiguity due to the presence of secondary 
lobes. 
3. Determination of an Acoustic Frontal Boundary 

To model the presence of an acoustic front two 
different sound speed profiles are introduced, one 50 km 
from the other. The parabolic equation model was utilized 
in this range-dependent problem with the position of the 
front (1.е., the range at which the second SSP is 
encountered) allowed to vary by +/- 20 km about the true 
position (step 1 km). Figure 3.10 provides an illustration 
of the SSP setup. The Bucker and the Heitmeyer methods 
correctly solve this problem with relatively narrow main 
lobes, as shown in Figures 3.11 and 3.12. However, the 
center of gravity method does not perform as well (Figure 
ЕКТІЗ) . Although it is able to locate the correct value, 
Significant sidelobes are present which could lead to an 


ambiguity if too small a detection threshold were chosen. 
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This last procedure is obviously not suitable for these 
conditions. 
4. Comments on the Estimators 
In order to completely appreciate the consistency of 
each of the previous estimators for the case where only one 
factor is unknown, it is important to study their response 
to a variety of configurations. 
a. Influence of the Number of Hydrophones 
In the problem of determining the depth of the 
SOFAR axis, the array was composed of 20 hydrophones. It is 
possible to run the same simulation by using only a fraction 
of the receivers. Plots of the Bucker estimator for four 
different numbers of hydrophones are shown in Figure 3.14. 
Although the difference is small when the number of 
hydrophones is reduced from 20 to five, the output of an 
array composed only of two receivers changes drastically. 
When the number of hydrophones is this small, the side lobes 
may have an amplitude of the same order as the main lobe, a 
Situation which leads to full ambiguity. The number of 
hydrophones is thus an important parameter which must always 
be more than some minimum value. This value is 
unfortunately dependent upon the actual problem and the 
depth of the array relative to the axial depth. Receivers 
which do not intercept much of the acoustic energy can be 
easily omitted but those which contain significant amplitude 


and phase information should be retained. 
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р. Influence of the Frequency 

Typically as a result of beamforming, the main 
lobe becomes narrower as the frequency of the array 
increases. The same phenomenon can be observed in Figure 
3.15, wherein the width of the estimator peak is also a 
function of the frequency. However, a trade-off exists 
between the desired resolution (width of the detection peak) 
and the computation time of the PE model which increases 
rapidly with frequency. Also if higher frequencies (kilo- 
Hertz range) were used, the signal would be limited Бу 
absorption which would result in lower signal-to-noise 
ratios. 

c. Importance of Array Position 

The depth of the array is of minimal importance 
when working in shallow water because the entire water 
column is nearly insonified. Such is not the case in deep 
water where shadow zones exist with relatively low signal 
levels. Based on the depth of the source, a first guess of 
the depth to position the array would be to place it where 
the signal may be expected to occur with a high level, for 
example, in the vicinity of the SOFAR axis or near a 
convergence zone. The choice will obviously be dependent on 


the profile shape. 
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C TWO PARAMETER PROBLEM (NOISY CONDITIONS) 
1. Localization of an Eddy 

The previous section has shown that the estimators 
are generally able to find the correct value in the case of 
a simple unknown and a noise-free medium. The same kind of 
simulation may be run when two parameters are to be 
determined. To test the ability of the various estimators 
to deal with a two dimensional problem, we will examine 
their ability to locate an eddy assumed to be present 
between a source and a receiving array. The sound speed 
profiles inside and outside the eddy are known. Thus the 
only unknowns are the borders of this perturbation of the 
sound field. An eddy 20 km in diameter is positioned 40 km 
to 60 km from the source. The replica fields are computed 
by scanning the limits from 35 km and 45 km for the border 
closest to the source, and from 55 km to 65 km for the 
farther boundary.  Replications at 1 km interval were made. 
We will thus try to match the simulated measured data with 
121 replica fields. Figure 3.16 presents the true location 
of the eddy in this simulation. 

Figure 3.17 shows the disposition of source and 
array With regard to the energy field for the true location 
of the eddy. The array lies on the SOFAR axis, almost 30 km 
beyond a convergence zone. From this plot, we can expect a 


high level of signal from the channel propagation and large 
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differences in phase due to multipath propagation (RR 
acoustic rays). 

Before processing, the measured data are imbedded in 
noise. This noise is introduced by the spatial autocorrela- 
tion matrix described earlier. The principal assumption of 
this simulation is that all correlation matrices are 
completely known. For an actual situation, this may not be 
true but it is still possible to estimate the total matrix 
of noise from the set of measured data. 

Because of its close similarity to the Bucker 
detection factor, the Heitmeyer function will be omitted 
from further analysis. The Bartlett and the Maximum 
Likelihood functions are introduced for these simulations, 
and it will later be seen that these two estimators are well 
Suited for conditions where the signal-to-noise ratio is 
low. 

2. Signal-to-Noise Ratio 

The signal-to-noise ratio 15 а parameter which 
depends on the relative powers of the signal and the noise 
at the array. Following the procedure of Ziomek (1985), the 
signal-to-noise ratio can be written in terms of the noise- 
free signal and the spatial autocorrelation matrix of the 


посе: 


ENR =e log), — (0482) 


The numerator of this expression may be interpreted 
as the summation of all the elements of the signal spatial 
correlation matrix. Similarly, the denominator represents 
the sum of the entries of the noise correlation matrix. 

The SNR was computed for several values of the noise 
power, of, and the correlation parameter, 8, that generate 
new values of KNi3 (када он ој The results аге 
presented in Table 3.1. The expression above shows that the 


SNR is reduced when the denominator of the argument 


NR NR 

increases, 1.е., when the double summation a КМ 1 - 15 
large due to a significant increase in Ес ен of the 
noise, =“, or a slow decay in the correlation between the 
hydrophones, :. By using the results presented in Figure 


3.1 and the analytic expression of the noise matrix given in 
Equation (3.1), the correlation of the noise will ре 
considered high when the factor © is less than 0.57 and low 


when it exceeds 2.0. 


TABLE 3.1 


SIGNAL TO NOISE RATIO AS A FUNCTION OF -? AND 


EAT 0.57 1.0 1 DN 22 
10719 13 ав 16 dB 20 dB 21 ав 21 ав 
5x 10-1? -1 dB 3 dB 6 dB 7 dB 7 dB 
10718 -7 ав -з ав о ав 1 ав 1 ав 
5 x 10-18 -21 dB  -17 dB -14 dB -13 dB -12 dB 
10717 -27 dB  -23 dB -20 dB -19 dB -19 dB 
10716 -47 dB  -43 dB -39 dB -39 dB -38 dB 
10718 -67 dB  -63 dB -60 dB -59 dB -58 ав 
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In the following sections, the performance of 
various estimators will be examined for several conditions, 
including different cases of noise power and noise 
correlation. 

3. Bucker Detection Factor 

A simulation using the Bucker detection factor was 
run for the case of eddy localization under both noise-free 
and noisy conditions. 

a.  Noise-free Conditions 

The simulation was run by setting the power of 
the noise, -?, to zero. The autocorrelation matrix of the 
noise is then just the null matrix and the total matrix 
reduces to one of signal only, as indicated by Equations 
ШИ!) and (2.29). A three dimensional plot and a contour 
plot of the detection factor are shown in Figure 3.18. The 
estimator is represented by a surface which has only a 
Single maximum positioned at the correct location of the 
eddy. This suggests that the Bucker method is able to 
determine the true location of the eddy in a noise-free 
environment. 

An interesting feature of the plot is the 
symmetry that exists around both diagonals of the contour. 


Moving on the principal diagonal, along the line: 


Ү = Х + 20 (3935 


29 


is equivalent to displacing an eddy with a constant diameter 
Of (26 сто The small intervals between the contour lines 
along this path indicate that the location of the eddy can 
be determined with reasonable accuracy, once we know its 
diameter. 

b. Case of Uncorrelated Noise 


Uncorrelated noise is generated by choosing a 


large value of ë. Im this case, the correlation falls off 
rapidly on either side of the noise matrix diagonal. For a 
large enough :, the noise field at one hydrophone is 


completely dissimilar to that at another hydrophone and the 
noise matrix becomes diagonal. Simulations were run with 
-= 10 and different values of 22. All yielded the same 
results as in Figure 3.19, which is seen to be identical to 
Figure 3.18. This similarity may be explained by recalling 


the definition of the Bucker detection factor when noises 





present: 
NR NR 
) ) TIED МА 
QUE c ke Ss JK 
BUCK = j-l k-17j (2.4) 
WO NM lJ 1/2 
() T R R ) ) КТ. Кт.) 
j=1 k=1⁄j 5 J 1=1 к=јиј 
where the total correlation matrix is given by 
KT- k = Коју + КМ- х (2.299 
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Or; 
КТік - К5ік % 22 ехр(-213-К|) (3.4) 


As can be seen, this expression only uses the 
cross terms of the matrix. When 8 is large enough, the 
second term in the right hand side of Equation (2.29) is 
almost negligible for every value of the noise power, 22%; 
thus the cross terms of the total matrix KT reduce to the 


cross terms of the correlation matrix of the signal KS. 
KTjk = К5-к, n (tor TES 


By changing the value of -?, the power of the 
noise is modified, but the new diagonal terms do not play a 
role in the calculations. The Bucker detection factor is 
thus insensitive to perfectly uncorrelated noise; in this 
case, the performance is exactly identical to the one in 
noise-free conditions. 

C. Case of Uncorrelated Noise 

Any combination of noise power, 22, and noise 
correlation, ;, yields a different pattern of the detection 
actor. In cases for which the spatial correlation of the 
noise is high, the Bucker method may still be able to 
maximize the detection factor at the correct location, but 


the absolute value of the peak will decrease as the 


S 


ambiguity surface becomes flatter. Figure 3.20 illustrates 
this type of behavior. 


Whenever both 2 апа 22 generate a low signal-to- 


noise ratio (large power, c?, or small correlation 
parameter, 8), the procedure fails and the localization of 
the eddy becomes impossible (see Figure 3.21). We will 


quantify the effects of c? and 8 on localization below. 
4. Bartlett Estimator 

The same simulations as above were run using the 
Bartlett estimator under noise-free and noisy conditions. 

a.  Noise-free Medium 

In a generic noise-free environment, the 

performance of the Bartlett estimator is similar to the 
Bucker detection factor, as shown in Figure 3.22. AS 


expected, the same symmetry along the diagonals is still 


present. 
b. Correlated Noise 
The performance of the Bartlett estimator 
changes significantly аз сб and vary. Figure 3.28 
provides an example of the plot for 22 = 10715 апа ё = 1. 
where the true location is found. Figures 3.24 and 3755 


illustrate failures of the method due to a weak signal-to- 
noise ratio brought about by strong noise and highly 
correlated noise, respectively. It will later be 
established that the usual characteristics of the noise in a 


deep ocean do not generally lead to this kind of ambiguity. 


32 


5. Maximum Likelihood Method 
The MLM estimator was calculated using the same 
noise conditions as for the Bartlett function. An 
examination of noise-free conditions is not possible because 
ef the singularity of the correlation matrix. Equation 
(2.31) shows that the expression of the Maximum Likelihood 
estimator requires the calculation of the inverse matrix 


KT71, 
MIM = (и: KT71 wy71 (01927 


When noise 15 absent, the matrix KT reduces to the 
correlation matrix of the signal KS which is generally 
singular. 
a. Slightly Correlated Noise 
When -? = 10719 апа — 1.7, the method gives 
better results than the Bartlett estimator with relatively 
low side lobes (compare Figure 3.26 with Figure 3.23). 
b. Strongly Correlated Noise 
With a more highly correlated ambient noise 
(> = 0.17), as depicted in Figure 3.27, it is still possible 
to obtain a correct location of the eddy. By comparing this 
plot with Figure 3.25, we note that the Bartlett estimator 
was unsuccessful in this case. Nevertheless, even for the 
MLM technique, a very low signal-to-noise ratio will result 


in a failure as shown in Figure 3.28. 
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IV. ANALYSIS OF THE SIMULATION 


From the previous simulations, one sees that the 
performance of each of the estimators varies significantly 
under varying signal-to-noise ratios or  source/receiver 
geometries or sound speed variations. For example, the 
center of gravity method was shown to be the best in 
locating the SOFAR axis and determining the surface sound 
speed. In contrast, this procedure was the least successful 
in the acoustic front localization problem. Therefore the 
efficiency of an estimator does not depend only on the type 
of ambient noise but also on the particular problem being 
solved. In order to continue focusing on a realistic 
problem, the eddy localization problem will be studied in 
greater detail. The condition of mismatching will be 


treated separately. 


Ax COMPARATIVE STUDY OF THE ESTIMATORS 
1. СЕ 
To facilitate comparison of the performance of the 
different methods, we will calculate the joint central 
moments of the different estimators. The joint «central 
moment provides information on the spread of the detection 
factor about the mean. The smaller the moment, the better 


the estimate of the parameter of the ocean. 


Gz 





Using the definition of the central moment of a 


multiple random variable as defined by Peebles (1987): 


илк = / /(х-Х)П (у-ү)К f ах а 2 
"nk / (х-Х) (Y-Y) ху (Х,У) х ау (4.1) 
where: 
Х, Y = means of the random variables X and MN 
fxy - joint probability density function of X and Y, 
п, К = orders of the central moment. 


The spreading factor for the MLM method is defined 


as the second order central moment of the function MLM(x,y): 


S= , MLM(X,y) (х-40)2 (у-60) 2 ах ду (485) 
where x and y are the boundaries of the eddy we are looking 
Гот. 

Since we are dealing with a discrete search among 


parameter values, this pseudo variance has the form: 


NR NR 


с 


S- ) } MLMj; (44+i-40)2 (54+j-60)* (ae) 
ізі 1-і 
where NF denotes the number of possible values of i and 3. 
Equation (4.3) indicates how the estimator spreads 
around the true frontal boundary values x = 40 and y = 60. 


It is a global measure of the estimator performance, not 
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just a measure of the main lobe width. A large value is 
associated with significant spreading which indicates a poor 
performance, even if the true position is actually found. 
The value of the spreading factor will be also large when 
the ambiguity surface has a large mean value and a small 
amplitude (case of significant noise).  Analogous spreading 
factors may be defined for the Bucker and the Bartlett 
functions. 
2. Efficiency of the Estimators 

The spreading factor defined above has been computed 
for the different combinations of -? and - shown in Table 
3.1. The results are presented in Table 4.1 for the Bucker, 
Bartlett and MLM methods. From this table it is possible to 
compare the efficiency of each technique in a variety of 
environments. The values of the correlation factor : have 
been chosen to stay consistent with deep water measures. 
The range of noise powers, s produced signal-to-noise 
ratios between -67 dB and +21 dB. 

In order to be consistent in comparing the spreading 
factor S of the three schemes, it is convenient to modify 
the expression for the Bucker detection factor as defined by 
Equation (2.4) by dividing it by its largest value among all 
the replicas. The maximum of this new function will always 
be one in all situations of noise, as the Bartlett and the 


MLM estimators. 
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TABLE 4.1 


SPREADING FACTOR S OF THE BUCKER (NORMALIZED), BARTLETT 
AND MLM ESTIMATORS AS A FUNCTION OF c? AND: 
E 0.57 GIO ДО 2.00 дао 
SNR=+13dB SNR=+16dB SNR=+20dB SNR=+21dB SNR=+21dB 
10719 19,000 18,973 18,950 18,944 8” ӘЛІ 
25,924 25 902 25,884 25,880 25,877 
114 139 142 140 139 
SNR--1dB  SNR=+3dB  SNR=+6dB  SNR=+7dB  SNR=+7dB 
5x10719 19,280 E Moos 19,004 19,090 
26 471 26,364 26,272 26,250 26,239 
567 654 708 698 691 
SNR--7dB  SNR--3dB  SNR-OdB SNR--1dB  SNR=+1dB 
10718 19,632 19,364 Ше 033 19,078 19,050 
27,119 26,935 26,754 "Tribe 26,688 
ms 1582 еа ые с 1. 377 
SNR=-21dB SNR=-17dB SNR=-14dB SNR=-13DB SNR=-12dB 
5х10718 22,400 Di GS 19,946 Ое а 19,532 
32,328 212300 30,463 30,258 30,153 
5,499 6,657 6,766 6,670 6,603 
SNR=-27dB SNR=-23dB SNR=-20dB SNR=-19dB SNR=-19dB 
mee! 625,773 23,188 20,950 20,407 Oo 
38,256 36,371 ЕЕ 34,386 34,190 
10,630 207/312 ВЕ 12,687 12,570 
SNR=-47dB SNR=-43dB SNR=-39dB SNR=-39dB SNR=-38dB 
io : 54,747 BL зала 30,280 
s 85,126 78,315 76,655 75,805 
66,998 yo 147 70,060 69,050 68,404 
, ,SNR--67dB SNR--63dB SNR=-60dB SNR=-59dB SNR=-58dB 
10 re = ee ex s 
со ° 4 B 124,781 
D oc oc 9 124,584 
Note 1: The signal-to-noise ratio is indicated for each 


entry of the table, followed 
represent, respectively, the 
and MLM spreading factors. 


by a triplet of numbers which 
Bucker (normalized), Bartlett 


Note 2: The value means that the estimator is unable to 
detect the true location of the eddy. 
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Table 4.1 provides a good illustration of mie 
performance of the estimators under several conditions of 
noise power and correlation. Figures 4.1, 4.2 and 4.3 show 
logarithmic plots of the spreading factor S versus the noise 
power 02; each curve represents a value of the correlation 
parameter 3. It is thus possible to determine the value of 
S for any pair of co^ and 8. As is obvious from Table 4.1, 
the comparative performance of each method is mostly a 
function of the signal-to-noise ratio of the measure. 

For SNR less than -50 dB all methods fail to detect 
the true location of the eddy. Nevertheless, we observe the 
case 52 = 10715 апа :- 2.2, where the Bartlett and MLM 
estimators indicate two different maxima at one. Even in 
this case, the amplitude of the functions is so small that 
the spreading factor S is very large. 

When the SNR is about -40 dB, the Bucker method is 
the most efficient estimator, albeit a weak one, as the 
Spreading remains significant. The superiority of the 
scheme increases moreover when the correlation of the noise 
decreases. 

For other SNR and correlation values, the MLM method 
is generally the most efficient. When the SNR = O GB or 
greater, the advantage of the MLM scheme iS obvious with 
respect to the other two functions. One also notes that the 
MLM estimator is well adapted to resolving highly correlated 


noise situations; in such cases, simulations show that the 
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ШЕЛ ОБ the main lobe becomes narrower but that the mean 


component of the surface increases. 


В. MISMATCHING CASE 
Mismatching occurs when a parameter used in simulation 

of the replica fields has been incorrectly estimated and is 
different from the one that created the true data. In the 
eddy localization problem, this defect may be introduced in 
several ways: 

- a wrong measure of the source frequency, 

- inaccurate estimation of the source or array position, 

- inaccurate estimation of the source beamwidth, 


- insufficient knowledge of the bottom 1055 
characterization, 


- oversimplification of the SSP. 

This list is not exhaustive and one must keep in mind 
that perfect matching almost never exists due to the 
impossibility of any acoustic model to solve the true 
acoustic wave equation in the real ocean. Because 
mismatching prevents a close likelihood between the actual 
and replica fields. It can be thought of as an additional 
noise which the correlation matrix has not taken into 
account. Mismatching thus results in a degraded estimation 
of the total spatial correlation matrix KT. The next 
section studies in greater detail two cases where 
mismatching is created by a change in the bottom loss 


parameterization and where the borders of the eddy are 
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smoother than what we have used previously to simulate the 
real measured data. 
1. Change in Bottom Loss Parameterization 

All simulations were conducted with the assumption 
that the bottom was fully absorbing. We will now consider 
the bottom to be a perfectly rigid surface with total 
reflection and examine how this new treatment modifies the 
llkelihood functions. 

In the absence of noise the difference in phase at 
each receiver of the array for both the perfectly reflecting 
and fully absorbing bottom conditions is depicted in Figure 
4. 4. When the bottom 15 treated as a perfect reflector, 
perfect matching will not occur because all the replica 
fields are constructed based on the full absorption 
assumption. 

Figures 4.5 and 4.6 show the result of a simulation 
using the MLM estimator with characteristics -2 - 10717 апа 
5 = 1.7 to represent the situations of no mismatching and 
mismatching, respectively. When mismatching occurs, the 
amplitude of the peak decreases due to an increase in the 
mean value of the likelihood function. The secondary lobes 
also become larger. It is important to note that the 
degradation observed in Figure 4.6 does not imply that 
treating the bottom as a perfect reflector is less correct; 


it only means that the matching was done improperly; one 
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must therefore be consistent in modeling the environmental 
input variables. 
2. Oversimplification of the Actual Eddy 

The previously described eddy localization 
Simulations were run with an excessive simplification of the 
true medium. The simulations assumed that there were no 
horizontal gradients of sound speed in and outside the eddy 
boundaries, 1.е., the change of SSP occurred almost 
instantly at the borders of the eddy. In an attempt to be 
more realistic, the next simulation was conducted after 
adding four intermediate SSPs between the two previously 
utilized profiles (Figure 4.7). The first intermediate SSP 
was introduced 4 km before the border of the eddy and the 
next ones added every kilometer thereafter. Actual signal 
values were generated with this smoothed baroclinicity, and 
replica fields were computed as before with the simplistic 
three-profile model. Noise was omitted from this simulation 
in order to better appreciate the effect of this 
mismatching. Results for the Bucker method are presented in 
Figure 4.8. Comparing this plot with Figure 3.18 we see 
that the main peak decreases but localization of the eddy 
remains possible. The implication of this simulation is 
that identification (location) of strong frontal boundaries, 
such as the north wall of the Gulf Stream or ice edge 
fronts, could be fairly exact but weaker, open ocean 


mesoscale eddies may pose more of a difficult problem 
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(assuming the same number of profiles are used to estimate 
the replica fields). 

The above examples suggest that in cases where the 
signal-to-noise ratio is quite low, it is розз1 Бей 
mismatching to hide the true location of the maximum and so 
possibly lead to a failure of the procedure. As Os i 
possible, mismatching must be avoided by a comprehensive 
knowledge of the parameters used in the replica fields 


calculations. 


CPC COMMENTS ON THE PROCEDURE 
As we were only interested in using matched field 
processing in acoustic tomography, many simplifications have 
been made to run the simulations. Although the procedure 
seems to be applicable and efficient in most cases, it is 
necessary to test its applicability to more complex 
problems. 
1. Environment 
It was assumed in all the simulations that only one 
or two parameters were unknown, for example, the surface 
sound speed or the eddy location. For actual oceanic 
Situations, many more properties can be expected to vary in 
Space and time. Extending this study to a heterogeneous 
medium is theoretically feasible but vastly increase ШЕ 
number of unknowns. 
Modeling a shallow water environment has not been 


considered as а possible mechanism to speed up the 
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calculations even though matched field processing remains 
possible in this kind of environment. Several limitations 
are apparent. A better Knowledge of the bottom structure is 
required. The ambient noise is moreover quite complex in 
coastal waters and its spatial correlation matrix would be 
difficult to model. As bottom interaction is important for 
estimator calculations (see the mismatching case), it is 
possible to consider the bottom loss as an unknown parameter 
and attempt to determine it through matched field 
processing. To examine this special case, the replica 
fields are generated using nine different bottom loss curves 
and one attempts to find the actual bottom loss curve. 
Figure 4.9 illustrates the different bottom loss curves that 
have been used to create the replica fields in shallow water 
(300 m). The maximum of the likelihood estimator occurs 
when the replica bottom corresponds to the actual loss. 
However, since the various curves are so similar in shape 
and because the bottom loss has only a weak to moderate 
effect on the transmission loss, the correct bottom 1055 
curve is not sharply defined. This implies that a correct 
specification of bottom loss for a low loss bottom is not 
required; not so for a high loss bottom. 
2. Model Consideration 
a. Resolution 
In the typical target location problem, the 


medium parameters are generally considered constant. It is 


П 


only necessary to run the PE model once to compute the 
pressure field at different distances. However, when the 
SSP becomes the unknown in the tomography problem, the PE 
model must be run for each replica. If we wish to find the 
correct shape of the SSP from 0 to 300 m in the seasonal 
thermocline with a resolution of 1m, the model needs to be 
run 10200 times if the sound speed at each depth can have 
ten different values. This simple example illustrates the 
trade-off between resolution and computer time. In the 
problem of eddy localization, where the limits of the 
perturbation were allowed to vary over 11 values, the PE 
model was run 121 times. For a determination of more than 
two unknowns, the basic theory of the estimators is still 
valid but the representation of the ambiguity surfaces 
becomes unachievable due to limitations in computer run 
time. 
b. Noise Approximation 

The correlation matrix of the ambient noise was 
modeled by a symmetric matrix decaying exponentially around 
the diagonal. This approximation is rather poor, even in 
deep water, because the power of the noise is assumed to be 
the same at each hydrophone. A better simulation would have 
been to consider the matrix of a noise which is, like the 
Signal, a solution of the wave equation. Our matrix Ше 


symmetric even though the actual matrix of the complex noise 


та 


needs to be hermitian, because the cross-correlations are 
complex. 
Sew Correlation Matrix 

If the spatial correlation matrix of the noise were 
known, it could be introduced in the replica fields and we 
could deal with it as with a noise-free problem. 

In practice, it 156 not possible to separately 
compute the noise and the signal correlation matrices. The 
total matrix needs to be estimated from the noisy signal at 
the hydrophones. Several techniques are available, as for 
example the Fourier method recommended by Johnson (1982). 
As the matrix becomes only an estimate; we should expect a 
slight mismatching and then decreasing efficiency of the 
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Ve- CONCEUSTON 


Matched field processing has been shown to be an 
efficient way to solve the inverse problem in ocean acoustic 
tomography when the ocean can be characterized by a few 
parameters. The estimators which have been used (Bucker, 
Bartlett, Maximum Likelihood, etc.) were generally robust 
enough to find the actual sound speed field of the ocean 


under usual conditions of noise power and noise correlation. 


A. SUMMARY OF THE RESULTS 

For noise-free conditions or high signal-to-noise ratios 
(SNR), all the estimators are able to correctly determine 
the actual unknown parameter of the medium. The Bucker 
detection factor was shown to be the best function when the 
SNR was moderate. 

The efficiency of the various estimators, illustrated by 
their spreading about the true value, decreased in cases of 
low SNR introduced by a large power or a high spatial 
correlation of the noise. The Maximum Likelihood method was 
shown to be the best scheme when the ambient noise was 
highly correlated because its spreading was less sensitive 
to the degree of spatial correlation than the other 
estimators. 

The effect of mismatching, when introduced in the 


simulations, generated a decrease in the efficiency of the 
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methods. Analysis of several degrees of mismatching 
indicated that unambiguous results can be expected from 
matched field processing provided that the parameterization 
of the medium is exact enough to generate consistent replica 


pressure fields. 


B. WEAKNESS OF THE SIMULATION 

In order to deal with reasonable computer times, only 
the cases of one or two unknown parameters were studied. 
For the same reason, the actual acoustic pressure field was 
compared with only a few replica fields. This limitation 
leads to moderate resolution in the results which could 
easily be improved by the generation of more replica fields. 

The weakness in modeling the noise field has already 
been emphasized in Chapter III.A. Further simulations 
should be done with a noise correlation matrix, KN, that has 
actually been obtained from measurements at sea. Further 
work using a more realistic parameterization of the ocean 
should be done, e.g., with empirical orthogonal functions 


(EOF) and using the technique to estimate EOF coefficients. 
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APPENDIX 


FORTRAN 77 PROGRAM USED IN THE SIMULATION 





PROGRAN ESTINA 
ACOUSTIC TOMOGRAPHY USING NATCHED FIELD PROCESSING 


THIS PROGRAM DRAWS THREE KINDS OF DETECTION FACTOR IN 3D 
IT COMPARES BUCKER,BARTLETY AND NLM METHODS 


THE NOISE IS INIRODUCED BY ITS CORRELATION MATRIX 


THE NOISE HATRIX IS PROPORTIONAL TO IDENTITY HATRIX IN CASE OF 
UNCORRELATED NOISE BUT IS IN GENERAL HERHITIAN MATRIX 


WIE EXACT COR sel) PHASES AND MAGNITUDES ARE READ ІМ THE 
FILE EXACT DATA 
THE COMPUTED PARAHRTERS ARE READ ON TUE FILLE NEAR DATA 


ОШ ООП АШ СУЫ до МАХТИЦИ OF 20 RECEIVERS DUE TO TIF USE OF THE 
PARABOLIC EQUATION FRON THE FEARG PACKAGE 


ES IS IE COUPLEX MATRIX OF MEASURED РАКАПЕТЕКӘ 

A IS THE NATRIX OF GUESSED PARANETERS WE WANT TO MATCH 

NR IS THE NUMBER OF RECEIVERS (CHAXINUM ay. 

NF IS TIE COMNNON NUNBER OF POSSIBLE VALUES FOR 2 UNKNOWNS 
THE NUNBER OF FIELDS WE NATCH 1S ACTUALLY МЕНЕ 


REAL lPIIASE (20) , PHAGHI (20) XO). жар DF(11,11),DFCGONT(11,11) 
REAL BART 11,11 „ЕМЪМ(Т1 10, HCONT(11,11) ВРАТ 
COMPLEX A 20; 20) Hsia 20) CORC1I. 11), hP(20, 20) 
COMPLEX W(20),CSQM, DETÉRM. ` 
REAL KN(20, 20) , HORE 
99 DEFINES THE FILE OF NEASURED(EXACT) DATA 
98 DEFINES THE FILE OF DATA WE WANT ТО МАТСИ 
NR=20 
NF=11 
READ TI THE EASURED VALUES IN FILE EXACT DATA 

READ (95, 600) PHAGNI(I),PPHASE(I) 
CONTINU 
FORMAT( — E11.4,2X,F7.4) 
COMPUTE TIE ELEMENTS OF MATRIX KS 
DO 2 J=1,NR 

DO 3 K=1,NR 

KS (J K)-PHAGNI (J)* PHAGNI (K)*EXP(CHPLX(0. , PPIASE(J) - 
PHASE (K))) 

CONTINU 

CONTINUE 


222) 


..4 0 0 Q 52... 6 0 Ч 0 ө 6 есе е е p. U U AC RCE CC ге песо пе е 9.9.9 9. ә өө T а 


С DEFINE THE CORRELATION MATRIX KN OF THE CORRELATED NOISE 
PRINT*,' VALUE OF SIGHA2 ? 
READ* , $ [GHA2 | 
PRINT*,' VALUE OF BETA 7 
READ* , BETA 
DO 4 'J-1,NR 
DO 5 Î=1,NR 
KN(1, JJ=SIGNA2* EXP (-BETA*ABS(I-J)) 
CONTINUE 
CONTINUE 


сл 


9$: 9 9g. т те оф оле cR o4 ро соо ое Е 


COMPUTE THE TOTAL COVARIANCE NAVTRIX INCLUDING THE NOISE 
ци QUO ДР NOISE AND THE SIGNAL CORRELATION MATRICES 
О = 


DO 7 Î=1,NR 
KT(I,J)-KS(I, J)ENCI, J) 
CONTINUE 

CONTINUE 


rere 


ар Гер өр Гель) 


BEGINNING OF HAIN LOOP FOR THE NF RUNS 


DO 101 ¿INDE5EZS 1 55 
DO 100 JNDEX=1,NF 
C READ PARAMETERS OF EXPECTED FIELD IN FILE NEAR DATA 


DO 8 I-i,NR 

READ(98,600) PHAGNI(T),PPHASE(L) 
8 CONTINUE 
С COMPUTE ELEMENTS OF HATRIX A 

DO 9 J=1,NR 
DO 10 K=1,NR 
A(J,K)=PNAGNI (J)*PMAGNI(K)*EXP(CMPLX(0., PPHASE(J)- 
* РРЛАБЕ(К))) 

10 CONTINUE 
9 CONTINUE 
0660 СС ОС СО О ОСС БЕБІ 
Ç DUCKER DETECTION FACTOR C 


СССССЕССЕССССССССЕЕСЕСССЕС обе 

CDF(INDEX,JNDEX)=CHPLX(0.,0.) 
00 11 J=I.NR 
DO 12 R=1,NR 

PECK NG J) CDF ( INDEX 
к A(,K)*CONJG(KT(CJ, K)) 
CONTINUE 

CONTINUE 


СЕКЕ x IS ACTUALLY REAL DUE TO PROPERTY OF HATRICES 
A AND KI {RENCE WE KEEP ONLY THE REAL PART OF IT 
DFCINDEX,JNDEX)-REAL(CDF(INDEX, JNDEX) ) 


NORMALIZE THE DF FACTOR 
FACTOR=0. 
FACT=0. 
0 13 J=1,NR 
00 14 K=1,NR 
1ЕСК. МЕ. J)FAGTOREFACTORIACI K)SCONJG(ACI К)) 
IF(K.NE. J)FACT-FACT4KT(J, K)*CONJG(KT(J,K)) 
CONTINUE 
CONTINUE 


JNDEX)-CDF ( INDEX , JNDEX) 4 


=— М 


(0 QC 


— — 
(9% ем 
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MOI Re СОЈ 


CI ANON 
~ со 


Jm 
15 


FACTOR-SQRT(FACTOR) 
FACT=SQRT (FACT 

DF (INDEX , JNDEX Е 
DF(INDEX, JNDEX)-ARS (DF ( INDEX , JRDEX ) ) 


CONTINUE 
CONTINUE 


NORMALIZE BY THE LARGEST VALUE 
BHAX=0. 
DO 67 I-1,NF 
DO 68 J=1,NF 
IF(DF(L,J).GE.BNAX) BNAX=DF(1,J) 
CONTINUE 
CONTINUE 


DO 69 I=1,NF 
DO 70 J=1,NF 
DF(I,J)-DF(I,J)/BHAX 
CONTINUÉ 
CONTINUE 


DISSPLAY THE NATRIX OF NORMALIZED DETECTION FACTOR 
НАИВНА BUCKER FACTOR HAIRIX 


DO 15 І-І,МЕ 
WRITE(6,777) (DE(I,J) ,J- 1, NF) 
ҒОВМАТ(11(25,ғ4.23) 


CONTINUE 


C 
CCCCCCCCCCCCCCCCCCCCCCCCCC 


BARTLETT ESTLUALOR C 


С 
CCCCCCCCCCCCCCCCCCCCCCCCCC 


С 


an 


DIEA 


сл 


aoc 


— — 
зо 


REWIND 98 


BEGINNING OF NATN LOOP TOR THIF NF RUNS 
DO 201 INDEX-I,NF 
DO 200 JNDEX-I,NF 


READ PARAHETERS OF REPLICA FIELD 
DO 50 I-1,NR 

READ(98,600) FHAGNHI(I),PFHASE(1) 
CONTINUE 


DETERMINE THE NORMALIZED COMPLEX VECTOR W 
NORN=0. 
DO 16 I=1,NR 
МОВИЕНОВИ + РИАСНТ (Г) 22 
CONTINUE 
NORN=SQRI (NORH) 


DO 17 I=1 
кое 
CONTINUE 
CONPUTE THE BARTLETT FACTOR 
CSUN=(0. ,0. 
DO 18 1-l,N 
DO 19 J=1,NR 
CSUMSCSUM4CONJG (W CI) )*KT CT, J)*W(J) 
CONTINUE 


CONTINUE 
BART(CINDEX, JNDEX)-REAL(CSUIT) 


NR 
|. /NORH) *PHAGNI(CI)* ЕХРОСМРІХ (0. „РРИАБЕ (1) )) 
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200 
201 


CONTINUE 
CONTINUE 


NORHALIZATION BY TIE HIGHEST VALUE 


BARNAX=0. 
DO 20 I=1,NE 
DO 21 J= 


1 ,NF 
LE BART (1 ,J).GE. BARHAX) BARNAX=BART (I,J) 
CONTINUE 


CONTINUE 
DO 22 I=1,NE 
DO 23 J=1,NF 
BART(1,J)=BARI(1,J)/BARMAX 
CONTINUE 
CONTINUE 


DISSPLAY TUE MATRIX OF NORMALIZED DETECTION FACTOR 
PRINT, > BARTLETT FACTORS NAIRIX 
DS i E 
WRITE(6,777)(BART(S,J) ,J=1,NF) 
CONT I NUF. 


CCCCCCCCCCCCCCCCCCCCC 


C 


NLM ESTINATOR C 


ООСС СББ СРЕ РОБ ОРО 


рер 2 


6262 


ca C 


26 


©? 


с) м 


28 


300 
30 


REWIND 98 


IHVERT THE CORRELATION HATRIX KT 
САБ СИТИ НЕ KL NR ODETERID) 


BEGINNING OF МАТН LOOP FOR THE NF RUNS 
DO 301 JNDEX=1,NF 
DO 300 JNDEX-1,NF 


READ PARAMETERS OF REPLICA FIELD 
DO 51 I-1,NR 

READ(98,600) PHAGNI(I),PPHASE(I) 
CONTINUE 


DETERMINE THE NORIALTZED COHPLEX VECTOR W 


МОКМ=О, 

КАЈ Пи 
ПОКПЕНОКИ+РНАСНТ (Г) #2 

CONTINUE 

NORN=SQRT (NOR) 

DO 27 I= 


1,NR 
cohen -(l./NORH)*PHAGNI (1)* EXP(CIHPLX(0 S PphpHASpCID))) 


COMPUTE THE MLN FACTOR 
CSUN=(0..0. 
DO 28 I-],N 
DO 29 191 МВ 
JN=CSUN4+CONJG (W(1))*RTCL, JWI 
CONTE (#(1)) ( )*W(J) 
CONTINUE 
FHLHCINDEX, JNDEX)-REAL(1. /CSUM) 


CONTINUE 
CONTINUE 
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С NORNALIZATION BY THE HIGHEST VALUE 
FNLNAX=0. 
DO 30 1=1 МЕ 
DO 31 J=1,NF 
IF(FHLHCI,J).GE.FHLMAX) FHLHAX-FHLH(CI,J) 
au CONTINUE 
30 CONTINUE 


C 
DO 32 1= МЕ 
DO 33. J=1,NF 
ЕКТ, J)=FNLN(1,J) /FNLMAX 
33 CONTINUE 
32 CONTINUE 
С 
C DISSPLAY THE MATRIX OF NORMALIZED DETECTION FACTOR 


PRINT*, MLM FACTORS HATRIX' 
DO За 1=1,МЕ 
WRITE(6, ВН ЕТ, НЕ) 
34 ПРОДАЈЕ 
С 


D 
ша пепоспсссессссссоссесесесесессессссссессссексессссс 
С СОНРОТЕ THE SPREADING FACTOR OF THE ESTIMATORS C 
И 
WDF=0. 
WBART-O. 
WEHLH=0 . 
DO 650 I-1,NF 
DO 651 J=1,NF 
WDE=WDF HDF (I nz КТ ) 2 (541 - 60) **2 
WBART-WBART4 BART is 1]-40. )***2*(54-]-60)9*2 
WEHLH-WFHLHTEHLH | Thi hh 41-40. обод. И ут? 


651 CONTINUE 
650 CONTINUE 

C 

С 


PRINI*,' WDF- ',WDF 
PRINI*,' WBART= | ,WBART 
PRINT*,' WFMLH- ',WFHLH 


END 
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